!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine check_state (ro, mx, my, mz, bx, by, bz, en, x, y, z, nx, ny, nz)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz

   real(8), intent(in) :: x(nx), y(ny), z(nz)
   real(8), dimension(nx, ny, nz), intent(inout) :: ro, mx, my, mz, bx, by, bz, en

   real(8) :: pr, vx, vy, vz, v2, b2, c2, s2, cax, cay, caz, cfx, cfy, cfz

!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx

      vx = mx(i, j, k) / ro(i, j, k)
      vy = my(i, j, k) / ro(i, j, k)
      vz = mz(i, j, k) / ro(i, j, k)
      b2 = bx(i, j, k) * bx(i, j, k) + by(i, j, k) * by(i, j, k) + bz(i, j, k) * bz(i, j, k)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en(i, j, k) - v2 * ro(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1

      if (ro(i, j, k) .le. 0.0d0) then
         ro(i, j, k) = eps
         mx(i, j, k) = 0.0d0
         my(i, j, k) = 0.0d0
         mz(i, j, k) = 0.0d0
         en(i, j, k) = pr / gmm1 + b2 / 2.0d0
         if (ab_density .eq. 0) then
            write(*, '(a, 1x, a, es10.2, 1x, a, es10.2, 1x, a, es10.2)') ' Warning: abnormal density at:',                         &
               'x', x(i), 'y', y(j), 'z', z(k)
            write(*, '(a, 1x, a, i10, 1x, a, i10, 1x, a, i10)') '                     index is: ',                                 &
               'i', i, 'j', j, 'k', k
         endif
         ab_density = 1
      endif
      if (pr .le. 0.0d0) then
         pr = eps
         mx(i, j, k) = 0.0d0
         my(i, j, k) = 0.0d0
         mz(i, j, k) = 0.0d0
         en(i, j, k) = pr / gmm1 + b2 / 2.0d0
         if (ab_pressure .eq. 0) then
            write(*, '(a, 1x, a, es10.2, 1x, a, es10.2, 1x, a, es10.2)') ' Warning: abnormal pressure at:',                        &
               'x', x(i), 'y', y(j), 'z', z(k)
            write(*, '(a, 1x, a, i10, 1x, a, i10, 1x, a, i10)') '                      index is:',                                 &
               'i', i, 'j', j, 'k', k
         endif
         ab_pressure = 1
      endif

      max_te_ro = max (max_te_ro, (pr / ro(i, j, k)) ** (5.0d0 / 2.0d0) / ro(i, j, k))

      c2 = gm * pr / ro(i, j, k)
      s2 = c2 + b2 / ro(i, j, k)
      cax = abs (bx(i, j, k)) / sqrt (ro(i, j, k))
      cay = abs (by(i, j, k)) / sqrt (ro(i, j, k))
      caz = abs (bz(i, j, k)) / sqrt (ro(i, j, k))
      cfx = sqrt ((s2 + sqrt (s2 * s2 - 4.0d0 * c2 * cax * cax)) / 2.0d0)
      cfy = sqrt ((s2 + sqrt (s2 * s2 - 4.0d0 * c2 * cay * cay)) / 2.0d0)
      cfz = sqrt ((s2 + sqrt (s2 * s2 - 4.0d0 * c2 * caz * caz)) / 2.0d0)

      max_speed = maxval ((/ max_speed, cfx + abs (vx), cfy + abs (vy), cfz + abs (vz) /))

   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine check_state
